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Abstract 

A new treatment of the classical beam torsion boundary value problem is applied. Using the p ~ version 
finite element method with shape functions based on Legendre polynomials, torsion solutions for generic 
cross-sections comprised of isotropic materials are developed. Element shape functions for quadrilateral 
and triangular elements are discussed, and numerical examples are provided. 


Introduction 

There are only a few known exact solutions for the torsion of beams. Among the exact solutions are those 
for elliptical, rectangular, and triangular cross-sections. The reason so few solutions exist is a problem of 
elasticity itself; analytical solutions to two- and three-dimensional boundary value problems for irregular 
areas or volumes are simply difficult to discover. Thus, approximate solutions have to be found. In this 
context, we extend the highly accurate ^finite element method to the two-dimensional boundary value 
problem of beam torsion. 

Based on Legendre polynomials, the p - finite element method offers exceptional convergence compared 
to the traditional h-version of the finite element method (Babuska [1]). In the p - version finite element 
method, the error in the solution is controlled by the polynomial order p : whereas in the h-ve rsion, 
the error is a function of the diameter h of the largest element. The advantages of p-FEM have been 
exploited in several areas, including elasticity, heat transfer (Smith [5]), and fluid dynamics. This is the 
first implementation of the method to the torsion problem. Numerical examples will show that its use in 
torsion problems is indeed valuable. The convergence of the solutions and its derivatives over the range 
from linear p = 1 elements to eighth-order p = 8 elements demonstrates the effectiveness of p-FEM to 
torsion. 


Classical Theory of Torsion 

In this section, we present the classical theory for torsion of beams for isotropic materials. If we consider 
a three-dimensional beam of length L with the cross-section shown in Figure 1 in the (xi, x 2 ) plane, at the 
end x 3 = L, a torque is applied. At the end z 3 = 0, the beam is constrained against rotation (translation 
of uj and u 2 displacements in the plane). Between the ends, the boundary T is free of stress. 

The state of stress in cross-section domain Cl must satisfy the equations of equilibrium 

divT = 0, (!) 

where for a linearly elastic, isotropic material, the stress tensor T is expressed in terms of the Lame 
constants A and p and the linear strain tensor E as 

T = A (trE) I + 2/xE. (2) 

The linear strain tensor E is given by 

E = i [gradu 4- (gradu) T ] , (3) 



Figure 1. Torsion problem domain. 
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where u is the displacement vector. In addition, the stresses must satisfy the condition on T of 


Tn = 0, (4) 

where n is the normal to the surface at any point on T. 

To solve the set of equations above, we assume a solution (Fung [4], p. 164) 

u\ ~ — ax 2 x 3 , u 2 ~ ax\X3, and 113 = a<p(x i,x 2 ), ( 5 ) 

where a is the angle of twist per unit length of the beam, and <p{xi, x 2 ) is the warping function. In (5), 
it is assumed that a is small such that the strains are strictly linear. Inserting (5) into (2), we find that 
the only nonzero stresses are 


Tl3 = \d^ ' X2 J and T23 = Ma \aZ * 

It should be noted that equations (6) identically satisfy (1) if we introduce the relations 


Ti3 ~ ftOt ~ 


and T 23 = 


where ip{x\, x 2 ) is the stress function. 

Differentiating the first equation in (6) with respect to x 2 and the second equation in (6) with respect 
to x\ and using (7), we have 

st is ( 8V A a 2 * 

ai7 = VaJIaJ7 “ ') = m 


d 2 <p 

dx\dx 2 


V = 


Subtracting equation (9) from (8) yields 


d 2 ip dhp 

dx\ dx\ 


which is the governing differential equation for the torsion problem for a linearly elastic, isotropic material. 
Recalling equation (4), we must meet the boundary condition, which is written as 

( dip dip \ /. .x 

73171! -f 7327X2 = ~ ) = 0 ‘ O 1 ) 


From Figure 1, we see that 


such that (11) can be rewritten as 


dx i 

and n 2 = — , 

d$ 


dip r 

— = 0 on r. 
ds 


Therefore, ip is constant on the boundary T, and we can set ip = 0 on T without loss of generality. 
The magnitude of the shear stress is given by 


T-yfa + n,.™ £ + Z 


anv where in the cross-section. In addition, the moment of the external forces at the end of the beam is 


= f {x\T 2 3 - x 2 T\3) dQ = — fia f 
Jn. Jci 


dip dip \ 

- + x 2 - — dQ ~ Da, 
dxi ox 2 J 
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Figure 2. Standard quadrilateral element geometry. 


where D is the torsional rigidity of the cross-section, or 

D= ~ ti L ( xi S +x2 £) = 2 »Jj dn - (ie) 

where Green’s theorem and ip = 0 on T have been used. 

In summary, the torsion problem is governed by the equations 

\7 2 ip = —2 in fi, il> = 0 on T, (17) 

with the stress magnitude given by (14) and the torsional rigidity by (16). 


Finite Element Solution of the Torsion Problem 

To solve equation (17), we apply Galerkin’s method. This method requires us to write the error residual 

R = V 2 ^> 4- 2 (18) 

and integrate it against a set of trial functions Ni over such that 

f RNidQ = f (V 2 i/> + 2) NidQ = 0. (19) 

Jn Jn 

The first term in equation (19) can be expanded through integration by parts as 

i vl * N ‘ dn = -/„(£;£: + Sir ) dn + L N, ^i n ' iT + L Ni %k n ' dT: (20) 

where the boundary terms vanish if (13) is used. 

If we write the stress function ip as 

'l’ = Y,a j N j , ( 21 ) 

>=1 

where aj are constant coefficients, then (19) can be rewritten as 

/ + = 1 2NidQ - {22) 
Jn \ 9 x\ dxi dx 2 0X2 ) Jn 


In terms of a system of equations, we have 


[K}{a} = {F}. 


The parameter N 3 is the number of trial shape functions Nj (or Ni) to be used in the solution. 
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Shape Functions for Quadrilateral Elements 

Before actually stating the shape functions used for the two-dimensional quadrilateral element, we should 
examine the geometry used. Figure 2 shows the quadrilateral element to be examined. The element has 
4 corners, 4 edges, and one interior, each of which must have shape functions associated with it. The 
corners in Figure 2 are labeled Cl to C4, and the coordinates of the corners in a (C, v) coordinate system 
are located at locations (±1,±1). The origin (C,» 7 ) = (0,0) is located at the center of the quadrilateral 


element in the parametric space. 

There is one mode associated with each corner for polynomial levels of p = 1 

and up. These shape 

functions are given by the functions 

O' 

II 

£ 

(24) 

where the C in (24) denotes a corner 

mode, and i refers to the corner number. 

Explicitly, the corner 

modes are given by 

Nci = ^(1 “ 0(1 “ v)> 

(25) 


Nci — ^(1 + C)(l “ 

(26) 


Nc3 = 4(1 + 0(1 “I” 

(27) 

and 

Nc4 = 4(1 - 0(1 + *7)- 

(28) 


Starting with p = 2, edge modes for the quadrilateral element are prescribed. For each additional 
polynomial level, four modes are added to the total number of degrees of freedom for the element, 
corresponding to the number of edges on the element. Thus, for any given polynomial order, we have 
4 (p— 1 ) edge modes, defined by 

N^v) = mv), (29) 

where E refers to the fact that this is an edge mode, p is the polynomial order of the element, and i is 
the edge number. These modes are written in terms of integrals of the Legendre polynomials such that 

Njg = \{l-r,)MQ, ( 3 °) 

N Ei = 5(1 + 0<M»?)> ( 31 ) 

N$ = 7j(l + h)<MO, (32) 

and 

tfis? = 5(1-0*07), ( 33 ) 

where the functions 0 , are defined by 

M*) = f_y P '~^ dX ' (34) 

and the functions Pj(x) are the Legendre polynomials of order j. 

For polynomial levels of degree p = 4 and up, internal modes become evident in the quadrilateral 
element. The internal modes are restricted to nonzero values inside the element (not at the corners or 
edges) . The internal mode for a p level of 4 is given by 

n { 4 0) = (1 -C 2 )(1 -n 2 ), P>4 (35) 

where the the subscript refers to the polynomial level, and the superscript denotes that this is an internal 
mode. 
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p 

Corner Modes 

Edge Modes 

Internal Modes 

Total DOF 

1 

4 



4 

2 

4 

4 


8 

3 

4 

8 


12 1 

4 

4 

12 

1 

17 

5 

4 

16 

3 

23 

6 

4 

20 

6 

30 

7 

4 

24 

10 

38 

8 

4 

28 

i 15 

47 


Table 1. Degree of freedom chart for quadrilateral elements 



Figure 3. Standard triangular element geometry. 

For higher modes, we can write the shape functions in terms of the shape functions for p = 4. That 
is, the shape functions for p = 5 and up have kernels composed of the shape function (35). For p = 5, 
there are two additional shape functions, given by 

<^2 j = JVi 0) {Pi(O, Pi fa)}, P > 5 ( 36 ) 

where P*(a) are Legendre polynomials of order i. The notation used in equation (36) associates two 
different functions with two different shape functions. For example, in equation (36), we have two shape 
functions, N^} and < 2 3 written as functions of Pi(C) and Pi fa), respectively. 

For the higher polynomial levels, we can write the face shape functions like 

^6°{1,2,3} = <<2(0, Pzfa), Pl(C)Pl(h)}, P>6 (37) 

<{1,2, 3,4} = < ) {^3(C),P3(»7),P2(C)Pl(»l).Pl(C)P2(*7)}, P > 7 (38) 

anJ <{1.2, 3 . 4 , 5} = <<4(0, P 4 fa), P 3 (C)Pl(»?). A(C)ftfa), P 2 (C)P2(0}- P > 8 (30 

Now that the shape functions have been written, we can easily determine the number of equations 
that will be needed for the analysis. We shall refer to the number of equations associated with a solution 
variable as the number of degrees of freedom per variable. The number of equations associated with each 
polynomial order is given in Table 1 . 

Shape Functions for Triangular Elements 

The shape functions for triangular elements are slightly more difficult since the element contains a slanted 
edge, as shown in Figure 3. The element has 3 corners, 3 edges, and one interior, each of which must 
have’ shape functions associated with them. The corners in Figure 3 are labeled Cl to C3, and the 
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p 

Corner Modes 

Edge Modes 

Internal Modes 

Total DOF 

1 

3 



r 3 

2 

3 

3 


6 

3 

3 

6 

1 

10 

4 

3 

9 

3 

15 

5 

3 

12 

6 

21 

6 

3 

15 

10 

28 

7 

3 

18 

15 

36 

8 

3 

21 

21 | 

45 


Table 2. Degree of freedom chart for triangular elements 


coordinates of the corners in a (C, 77 ) coordinate system are located at locations (0, 0), (TO), and (0, 1) in 
the parametric space. 

There is one mode associated with each corner for polynomial levels of p = 1 and up, given by 


and 


Nci = 1 - C - i?. 
N C 2 = C, 

Nc3 = rj. 


(40) 

(41) 

(42) 


Starting with p = 2, three edge modes are added to the total number of degrees of freedom for the 
element. Thus, for any given polynomial order, we have 3(p — 1) edge modes, defined by 


tfiv = (!-<- v)MC), 
K ( Ei = nMO, 


(43) 

(44) 


and 


-•?)*(>?)• ( 45 ) 
For polynomial levels of degree p — 3 and up, internal modes become evident in the triangular element. 
The internal mode for a p level of 3 is given by 

ivf : ) = LiLiLa, P> 3 (46) 

where the functions L\ , L 2 , and L 3 are defined as 

Li = 1 - c “ I 2 = C, and L 3 = ( 47 ) 


For higher modes, we can write the shape functions in terms of the shape functions for p = 3. That 
is, the shape functions for p = 4 and up have kernels composed of the shape functions (46). For p = 4, 
there are two additional shape functions, given by 

^°{1,2} = ^3 0) { ^2 - L l)> *M 2L 3 - !)}• P > 4 (48) 

The remaining internal modes are given by 

12 3 , = .vl 0) {P 2 (I 2 - Li), P 2 (2T 3 - 1), Pi(Li - Li)Pi(2L 3 - 1)}. P > 5 (49) 

Ni°l 234) = N { 3 0) {P 3 (L 2 -L 1 ),P 3 (2L z -l) l P^L 2 -L 1 )P 1 (2L 3 -l),P 1 (L2-L 1 )P 2 (2L 3 -l)}, p>6 

(50) 
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Figure 4. Single-element meshes for full and quarter models. 
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Figure 5. Dimensionless torsional constant versus p-level for a square beam. 


1 V 7, { 1,2, 3, 4,5} 


Ni 0) {Pi{L 2 - L x ), P 4 (2L 3 - 1), Pz{L 2 - Ii)Pi(2I 3 “ 1), Pi(L 2 - Ii)P 3 (2I 3 - 1), 

P 2 (L 2 -L 1 )P 2 (2L 3 - 1)}, P>7 (51) 


and 

^8°{1,2,3,4,5,6} = ^ 0) {Ps(L 2 - Ii),P 5 (2L 3 - 1), Wa - Li)Pi( 2I 3 - IJ.Wa - U)P A {2L Z - 1), 

P 3 (L 2 -I 1 )P2(2L 3 -l), J P2(l2--Li)P3(2L3- 1)}. p > 8 (52) 

The number of equations associated with each polynomial order for triangular elements is given in 
Table 2. As can be seen from Table 2, there is not much difference between the total number of degrees 
of freedom associated with a triangular element and that of the quadrilateral element, shown in Table 1. 
Thus, it is preferable to use quadrilateral elements when defining geometry. 


Torsion of Beams with Square Cross-Sections 

In this section, we demonstrate the use of the p- version finite element method for torsion problems. If 
we consider a beam with a square cross-section, then we can write the exact solutions for the maximum 
torsional shear stress and torsional constant as (Sokolnikoff [6], p. 131-132) 


7~max — fXOtCL 


1-4 

1 


cosh7r/2 


^ (2n + l) 2 cosh(2n + 1 )tt/2 


7 
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Figure 6. Dimensionless torsional shear stress versus p-level for a square beam. 


and 


D = /ia 4 


1 64 ^ tanh(2n -f 1)7 t/2 

3 " ** \ (2 n + 1)5 

n =0 


(54) 


where a is the length of one of the sides of the cross-section, fi is the shear modulus of the beam material, 
and a is the angle of twist per unit length of the beam. 

To compare the analytical solutions (53) and (54) to the finite element solutions, we consider two 
different meshes. First, if we consider the left mesh of Figure 4, we have a single element representing 
the domain Q. Next, by exploiting the double symmetry of the cross-section, we can model one-quarter 
of the cross-section, yielding the second mesh in Figure 4. We note that in this mesh, on two sides we 
have dip/dn = 0. 

In Figure 5 we show the convergence of the normalized torsional constant D/pa 4 as the polynomial 
order of the solution is increased. The flat, horizontal line is the value given by equation (54). It should 
be noted that for the full model, there are no computed values of the normalized torsional constant until 
p = 4. The reason for this is that from p = 1 to p = 3, there are only corner and edge modes. Since 
ip = 0 on the boundary, then clearly the solution contains no active degrees of freedom, and the solution 
is forced to zero. However, at p = 4, we obtain the first internal mode, and its value is reflected in the 
computation of a non-zero torsional constant. At p = 4, the full model mesh only has the single active 
degree of freedom, yet its computed torsional constant value has less error than the quarter model for 
p — 3, which has more degrees of freedom. Thus, it appears that additional degrees of freedom are not 
terribly important in evaluating function values. 

However, when shear stresses, which are derivatives of the trial solution, are evaluated, the calculated 
flux values are somewhat dependent upon the order of the solution, as seen in Figure 6. Again, for 
p = 1 to p = 3, the one-element full model solution is identically zero since all of the degrees of freedom 
are constrained to zero from the tp = 0 restriction. Thus, the derivatives are also zero. However, with 
increasing p, we do see the existence of a derivative value for the full model, and over the entire p range, 
the quarter-model shows exceptional convergence. 


Torsion of an Elliptically Shaped Beam 

The simple cross-section of the previous section can be expanded to a beam with an elliptical cross-section. 
The ellipse provides yet another known theoretical solution against which the numerical finite element 
solution can be compared. In this case, the maximum shear stress due to torsion, which occurs closest to 
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Figure 7. Two element mesh for quarter-model of elliptical section. 
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Figure 8. Dimensionless torsional constant versus p-level for an elliptical beam. 


the major and minor axes of the ellipse, and the torsional rigidity are given by (Fung [4], p. 170) 

a 2 b 


— 2fia 


a 2 *4* b 2 


and 


D = \i 


7 ra 3 6 3 
a 2 + b 2 5 


(55) 


(56) 


where a and b are the major and minor axes of the ellipse. 

For the comparison, a simple two-element mesh representing one-quarter of the entire ellipse is em- 
ployed. One quadrilateral element and one triangular element are used, as shown in Figure 7. For 
an ellipse with a = 3 and b = 1.5, we find that for p > 2, there is good agreement with the normalized 
torsional rigidity D/fi, as shown in Figure 8. Likewise, there is rapid convergence of the normalized shear 
stress T max / jj, , as shown in Figure 9. The reason for the quick convergence is that the exact solution is 
given by (Fung [4], p. 169) 

*— (57) 

Since (57) is quadratic in x and y, the theoretical solution is quickly converged for p = 2, which is shown 
in the graphs for the torsional constant and the maximum shear stress. 


Torsion of a Generic Section with Cutouts 

Textbook solutions for the torsion problem can be demonstrated with finite element methods easily, but 
they are fairly uninteresting and show only that we can match solutions. To demonstrate the effectiveness 
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Element Order 

Figure 9. Dimensionless torsional shear stress versus p-level for an elliptical beam. 


v=0 



Figure 10. Finite element mesh of a generic cross-section. 

of the finite element method to torsion, especially with p-finite elements, we examine the torsion of the 
cross-section shown in Figure 10. The cross-section has symmetry about two axes, so only one-quarter of 
the cross-section is modeled. Along the outside of the section, the function ip is identically zero. Likewise, 
along the perimeter of the cutouts, the stress function is forced to zero. Along the straight edges, which 
represent lines of symmetry, we have the condition dip j dti — 0. The finite element mesh is a seven-element 
mesh using eight noded elements to achieve the curvature along the edges. 

The cross-section is excellent at demonstrating the torsion of multiply connected regions which would 
otherwise be left unsolved by analytical means. The torsion analysis was conducted with element shape 
functions from p = 2 to p = 8. The results for p = 1 are inconclusive since the p = 1 shape functions are 
corner modes, and the mesh requires that ip — 0 at the corners. The convergence of the torsional constant 
is shown in Figure 11, whereas the shear stress at a point in the left cutout is shown in Figure 12. Like 
the cases shown previously, we find rapid convergence for both the function value and the first derivative. 

Finally, an error estimate of the solution is plotted in Figure 13. Two separate mesh strategies 
were employed in this analysis. The first mesh is the same as shown in Figure 10. That is, the mesh 
remains the same, but the polynomial order of the solution is increased from p = 2 to p = 8. The second 
mesh strategy is to take the mesh of Figure 10 and subdivide the areas into smaller regions with uniform 
h-refinement. That is, each element is subdivided into n 2 elements, where n is the number of elements 
per side. This yields a uniform h-refinement and comparison curve to demonstrate the convergence of 
the new technique. The error in the mesh is plotted in Figure 13 versus the number of active degrees of 
freedom in the solution, and it is defined as 

e= \n mesh - Doc\, (58) 

where £) mejh is the torsional constant corresponding to the current mesh and D a 0 is the torsional constant 
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Figure 11. Torsional constant versus p-level of the solution. 



Figure 12. Torsional shear 


of the solution. 






NDOF 


Figure 13. Error estimates for the p-element and refined h-e lement meshes. 

computed for a mesh with p = 8 and a significant number of elements. In Figure 13, we note the 
exceptional performance of the p-version mesh. While the refined h mesh converges at virtually a constant 
rate, the p-element mesh, composed of just 6 elements, rapidly converges with the addition of higher order 
polynomials. 


Conclusions 

As demonstrated by the examples, implementing the p- version of the finite element method is beneficial 
in achieving highly accurate, converged solutions to the torsion boundary value problem. The flexibility 
of increasing the fidelity of the solution by changing the polynomial order rather then remeshing to find 
appropriate element size saves both time and effort without sacrificing numerical accuracy. Indeed, it was 
shown that for equal degrees of freedom, the p-method far outperforms the /i-method in the calculation 
of the torsional constant. While no results for the shear stresses were presented with the general cross- 
section, we would expect similar if not better characteristics, as the h- method is limited to at best first 
derivatives which vary linearly in the element space. However, for the ^method, the derivatives are of 
higher order, giving them the ability to match gradients better. 
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